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Abstract 

Thermostated tethered harmonic lattices provide good illustrations of the phase- 
space dimensionality loss AD which occurs in the strange-attractor distributions 
characterizing stationary nonequilibrium flows. We use time-reversible nonequilib- 
rium molecular dynamics, with two Nose-Hoover thermostats, one hot and one cold, 
to study a family of square heat-conducting systems. We find a phase-space dimen- 
sionality loss which can exceed the dimensionality associated with the two driving 
Nose-Hoover thermostats by as much as a factor of four. We also estimate the di- 
mensionality loss AZ>ft in the purely Hamiltonian part of phase space. By measuring 
the projection of the total dimensionality loss there we show that nearly all of the 
loss occurs in the Hamiltonian part. Thus this loss, which characterizes the extreme 
rarity of nonequilibrium states, persists in the large-system thermodynamic limit. 
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1 Introduction 



Aoki and Kusnezov stressed the usefulness of the "0 4 " model in studying the 
size-dependence of classical heat conductivity simulations [1]. They used this 
model to show that the entropy production associated with steady heat flow 
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becomes extensive, (oc A), as the system size A increases. When such heat- 
flow simulations are driven with two time-reversible Nose-Hoover thermostats 
[2,3,4,5,6,7] imposing two temperatures, "hot" and "cold", the steady-state 
phase-space distribution function becomes a multifractal object [8,9,10]. The 
information dimension of this multifractal object lies below the corresponding 
equilibrium one by an extensive (oc A) "dimensionality loss" AD. The dimen- 
sionality loss, like the closely-related thermodynamic entropy production S, is 
approximately quadratic in the deviation from equilibrium brought about by 
the temperature gradient VT: 

AD oc iV(VlnT) 2 . 

This quadratic dependence is the usual prediction of linear transport theory. 
Its evaluation quantifies the rarity of nonequilibrium stationary states. 

About ten years ago [11] we studied phase-space dimensionality loss in a fam- 
ily of two-dimensional shear flows. Those simulations were limited to relatively 
small systems. In that shear-flow work we were able to find systems for which 
the dimensionality loss barely, but significantly, exceeded the total dimen- 
sionality associated with the coordinates, momenta, and friction coefficients 
{q,p,C} °f the thermostated boundaries. 

At that time we stated that our data definitely showed that nonequilibrium 
stationary states can exhibit a reduced dimensionality, relative to equilib- 
rium, even when projected into the subspace occupied by purely-Hamiltonian 
degrees of freedom. This interpretation was controversial [12,13]. Many re- 
searchers had difficulty in accepting that nonequilibrium stationary states are 
typically associated with reduced dimensionality. In order to confirm that our 
interpretation was correct, we recently studied dimensionality reduction for 
stationary flows using the (Hamiltonian) 4 model. Because heat transfer us- 
ing this model requires only the simplest of boundaries — a single hot degree of 
freedom and a single cold one, in the simplest one-dimensional case — it turned 
out to be easy to get relatively large dimensionality losses [2]. 

A simple but convincing two-dimensional case involves the heat transfer from a 
single hot particle (with five phase-space coordinates {xH,UH,PxH,PyH,(H}) 
to a single cold one (with its own five coordinates) through a medium of 
A — 2 = 23 Hamiltonian particles (with 92 more phase-space coordinates, 
for a total of 4 A + 2 = 102 phase-space coordinates). Typical (x, y) particle 
trajectories in such a 25-particle system are shown in Fig. 1. In the case 
of A = 16 particles discussed in Ref. 2 the phase-space strange attractor 
turned out to have a dimensionality loss of 12.5, relative to the equilibrium 
distribution, exceeding the number of thermostated phase-space coordinates 
(5 hot and 5 cold) by 2.5. 

The present work has two goals, (i) characterizing the overall loss of phase- 
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Fig. 1. Geometry of a 25-particle two-dimensional steady-state system with a 
full-phase space dimensionality reduction AD = 21.6. The upper righthand par- 
ticle is "hot" and the lower lefthand one is "cold" . The remaining 23 particles obey 
conventional Hamiltonian mechanics. The trajectories shown were generated with 
hot and cold temperatures Th = 0.009 and Tc = 0.001 imposed by two Nose-Hoover 
thermostats with relaxation times of r = 5. 



space dimensionality AD, and (ii) estimating the corresponding loss AD^ 
in the Hamiltonian unthermostated part of phase space. This two-part study 
requires the numerical evaluation of the instantaneous Lyapunov exponents 
{A(t)}, as well as their associated vectors {5}. The required techniques are 
sketched in Sec. II. In Sec. Ill we describe the 4 model. In Sec. IV we under- 
take a systematic study of the number- dependence of the phase-space dimen- 
sionality loss AD. We consider a family of systems of the type illustrated in 
Fig. 1. The numerical work described in Sec. IV shows convincingly a \Jl/N 
deviation from the large-system limit, with an extrapolated dimensionality 
loss of as much as four times (40) the number of dimensions (10) associated 
with the thermal driving mechanism. Next, in Sec. V, we study the projection 
of this dimensionality loss into the purely Hamiltonian subspace. We find that 
nearly all of the loss occurs in that subspace despite the purely Hamiltonian 
form of the equations of motion there. We discuss the projection technique 
and the importance of rapid rotation in phase space [14] to an understand- 
ing of the strange-attractor's dimensionality reduction. The conclusions which 
follow from this work take up the final section. 
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2 Lyapunov Exponent Calculations 

The usual "Lyapunov exponents", {Xj}, are time averages, over a sufficiently 
long time for convergence, of "local" (instantaneous) exponents {Xj(t)}, 

A, = (A,(t)> . 

Here, and in what follows, we use angular brackets (...) to indicate long- 
time averages. Each of the instantaneous exponents has associated with it an 
"offset" vector 5j in phase space which describes the direction in which the 
growth or decay of phase-space separation is measured. If the offset vectors 
are infinitesimal in length, as in the present work, they can be replaced by a 
parallel set of unit vectors in "tangent space" . 

The time evolution of the offset vector directions is governed by a continuous 
Gram-Schmidt orthonormalization which forces the vectors to remain mutu- 
ally perpendicular and to evolve at fixed length. The two types of constraints, 
orthogonality and fixed length, can be imposed by a triangular array of La- 
grange multipliers, A^-, where 1 < j < i < n in an n-dimensional phase space 
[15]. When the usual Cartesian coordinates are used the orthonormal offset 
vetors rotate rapidly [14]. 

The usual Lyapunov exponents {A-,} are the long-time- averaged values of the 
diagonal Lagrange multipliers {Ajj} required to enforce the orthonormaliza- 
tion constraints: 

A, = (A,(t)> = (A,,) . 

At any instant of time the instantaneous Lyapunov exponents X(t) represent 
the orthogonal measurements of the n growth (or decay) rates determined by 
the dynamical matrix D. D, an n x n matrix, is itself a phase function. The 
ith row of D is made up of the n derivatives of the ith equation of motion 
with respect to the n phase variables. 

The eigenvalues and eigenvectors of the dynamical matrix are relatively com- 
plicated [16]. The eigenvalues are mainly rapidly varying complex-conjugate 
pairs, with the imaginary parts frequently vanishing, at phase-space singu- 
larities corresponding to parallel eigenvectors. Despite this complexity the 
orthonormal basis provided by the offset vectors makes it possible to measure 
smooth and well-defined growth rates. The time reversibility of the equations 
of motion guarantees related reversibility properties for D and its eigenvectors 
and eigenvalues. Nevertheless, the past-based Lyapunov vectors governed by 
D show a time-symmetry breaking intimately related to the second law of 
thermodynamics [3,8,9,11,14,16]. 

Consider a simple textbook [3] (Sec. 5.4) illustration of these ideas, the motion 
of a driven thermostated particle in one dimension. The motion takes place in 
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the three-dimensional (q,p,() phase space with three Lyapunov exponents: 



Ai 


= (Ai(0> = 


(An) 


A 2 


= <A 2 (f)> = 


(A22) 


A 3 


= (As(O) = 


(A33) 



The equations of motion (with unit mass, force constant, and relaxation time) 
are 

q = p; p = +l-(p; ( = p 2 - 1 , 
for which the long-time solution is an attractor: 

{(Z(*),P(*U(*)} — 0, +!,+!}■ 
The matrix of equation-of-motion derivatives D is 



^0 +1 " 



D = 



-c - P 

2j) 



On the attractor the momentum p and friction coefficient ( can be replaced 
by their limiting values, (+1,+1). The three tangent-space S vectors follow the 
ordinary differential equations: 



61 = D ■ Si - A n 5i 

62 = D ■ 5 2 — A 2 i5i — A 2 2#2 

S 3 = D • S 3 - A 31 Si - A 32 S 2 - A 33 S 3 , 

where the six Lagrange multipliers follow easily from the time derivatives of 
the six orthonormality conditions: 



2 r2 
" S 



S 2 3 = l; 



S\ ■ s 2 = s 2 ■ s 3 = S 3 ■ Si = . 

The steady-state time variation of the Lagrange multipliers, and the three S 
vectors are shown in Fig. 2. The three Lyapunov exponents, like their three 
instantaneous values {\j(t) = Ajj}, sum to —1: 



| + | + |, Al(t) + A2(t) + As(t) 



-1 . 



Note that the largest (time-averaged) Lyapunov exponent is in this case, 
corresponding to a phase-space displacement in the direction of the motion 
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Fig. 2. Bottom: Time dependence of the five nonvanishing Lagrange multipliers for 
a field-driven thermostated particle in one dimension, as discussed in Sec. II. The 
time dependence of the vectors {82, 83} for this same problem is shown at the top. 
The period of the oscillation is 2-^-^/4/7 though the phase has no significance in the 
stationary state. 

5i = (1,0,0). The remaining five Lagrange multipliers, as well as the vectors 
82 and <5 3 , oscillate periodically, with a period of 27ry / 4/7. 

An analogous computation for an n-dimensional set of vectors involves com- 
putational work of order n 3 . Nonequilibrium steady states generate chaotic 
attractors rather than the simple fixed point of the damped oscillator ex- 
ample. Nevertheless, the basic steps are the same: (i) propagating a phase- 
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space reference trajectory; (ii) simultaneously propagating offset vectors in 
the neighborhood of the reference, with Lagrange multipliers and/or Gram- 
Schmidt orthonormalization imposing orthonormality; and (hi) averaging the 
diagonal multipliers to find the Lyapunov spectrum. Let us now consider the 
many-body model for which we carry out such investigations. 



3 4 Heat Flow in Two Dimensions 

In the "0 4 " model we use here [1,2] we choose a square nearest-neighbor har- 
monic lattice with a quadratic Hooke's Law pair potential for neighboring 
particles % and j, 

0(r) = ^f(r — d) 2 ; r = \n — rj\ > . 

In addition, each particle is tethered to its lattice site with a quartic potential, 
^-Sr 4 . The quartic tethers have two nice consequences. First, they provide 
"external forces" , and so prevent momentum conservation and ballistic energy 
transport. Second, they can provide chaos, so that the dynamics can have one 
or more positive Lyapunov exponents. The two consequences together can 
give Fourier heat conduction, even in one or two space dimensions. The heat 
conductivity for the 4 model remains finite in the large-system limit, unlike 
many models, for which the conductivity vanishes or diverges in the large- 
system limit. 

Here we simulate nonequilibrium heat-conducting stationary states by impos- 
ing thermostating forces on two of the N particles. The thermostating forces 
are Nose- Hoover [3,4] feedback forces {— (p}, linear in both the time- reversible 
friction coefficients {(} and the momenta {p} = {(p X) p y )}. The thermostated 
equations of motion for the hot and cold particles (one of each) are: 

{mx = p x ;my=p y } 

{ Px = F x - CPx ; P y = F y - (p y } 

C(h or c) = [(p 2 /2mkT (H o r c))-l]A 2 , 

where r is the characteristic response time of the thermostat forces, {— (p}. 
The full phase space describing this iV-body two-dimensional system has 47V+ 
2 dimensions, with the extra two corresponding to the hot and cold friction 
coefficients ( H and Co 

For convenience in numerical work we choose the particle mass m, the spring 
constants k 2 and k 4 , and the nearest-neighbor lattice spacing d all equal to 
unity. The remaining parameters to set are the hot and cold temperatures, 
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which we arbitrarily choose equal to the values from Ref. 2, 

kT H = 0.009 ; kT c = 0.001 , 

and the thermostat relaxation times. We vary these times in the numerical 
work, but with the simplifying restriction that th and tq have a common value, 
which we denote as r. We use the classic fourth-order Runge-Kutta integrator 
throughout, with a timestep dt = 0.001 or 0.002. To avoid numerical errors we 
have generated and compared results from two fully-independent simulation 
codes, one written in Vienna and the other written in Livermore. 



4 Numerical Results — AD 

The numerical evaluation of the dimensionality loss AD is based on the con- 
nection between the Lyapunov spectrum and the dimensionality of the phase- 
space strange attractor. The Lyapunov exponents give the time-averaged rel- 
ative growth and decay rates of the principal axes of a comoving infinitesimal 
phase-space hypersphere (or "extension in phase"). Kaplan and Yorke con- 
jectured, evidently correctly [17], that a partial sum of these exponents (be- 
ginning with the largest one and proceeding toward the most negative one) 
changes from positive to negative when the (linearly-interpolated) number of 
terms in the sum is equal to the dimensionality of the phase-space strange 
attractor. This conjecture is "almost obvious". It is evident that the (hy- 
per)volume of any phase-space object with a positive sum of time-averaged 
Lyapunov exponents must diverge. Likewise a negative sum of time-averaged 
exponents indicates a vanishing hypervolume at long times. Any stationary 
process must generate an attractor which neither vanishes nor diverges. 

The main difficulty in computing Kaplan- Yorke information dimensions is the 
unfavorable time-dependence associated with constraining the phase-space 
offset vectors to remain perpendicular to one another. With Gram-Schmidt 
orthonormalization in an n-dimensional phase space n vectors, with n com- 
ponents each, must all be propagated in time for sufficiently long that the 
time-averaged growth rates have converged. The computational work in or- 
thogonalizing n(n — l)/2 pairs of n-dimensional vectors varies as the cube of 
the number of particles, so that present computer speeds and processor num- 
bers allow us to follow no more than a few hundred particles. In the present 
work we consider the simplest possible square systems of from 4(2 x 2) to 
144(12 x 12) particles. The relaxation time r is a free parameter. We choose 
it in the range 1 ... 8. Results become insensitive to r once r exceeds 6, which 
is comparable to the inverse Debye frequency of the lattice. Representative 
results for the largest Lyapunov exponent Ai and the dimensionality loss from 
the Kaplan- Yorke conjecture, AD, are given in Table I. 
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Table 1 



Representative data for the total phase-space dimensionality loss AD and the 
Hamiltonian projection AD^ as a function of system size and the Nose-Hoover 
relaxation time r. The largest Lyapunov exponent, Ai is also tabulated. The two 
boundary temperatures (imposed by a single hot and a single cold particle) are 
0.009 and 0.001 in all cases. The estimated error is ±0.4 for AD and ADn, and 
±0.0005 for Ai. 
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Our own previous work, on color conductivity and shear flow [11], strongly 
suggests deviations in the Lyapunov spectrum of order the inverse system 



width, yl/iV in two dimensions. The present results are roughly consistent 
with this finding though an even slower variation with N provides a compa- 
rable fit. We know of no previous systematic study of the variation of the 
Lyapunov spectrum with r. A cursory investigation shows that the dimen- 
sionality loss varies roughly as r _1 . Thus the dimensionality loss, with fixed 
boundary temperatures of 0.009 and 0.001, can be represented by 



The standard deviation affecting the last digits of the fit parameters are given 
in brackets. We compare the second expression with our numerical data in 
Fig. 3. The rather good fit indicates that the dimensionality loss persists in 
the large-system limit. We investigate the loss further in the following section. 




AD = 42.7(5) 



89(2) 13.2(5) 
VN t 
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Fig. 3. Comparison of the fit AD = 42.7 - (89/iV 1/2 ) - (13.2/r) to the simulation 
results for the dimensionality loss in the full 4N + 2-dimensional phase space. 

5 Hamiltonian Projection: AD n 

Ten years ago [14] we studied the rotation rate of the n phase-space offset 
vectors {Sj}. We found that the rotation rate increases very rapidly with sys- 
tem size, soon becoming very large relative to the Lyapunov exponents them- 
selves. This observation suggests that the instantaneous growth and decay 
rates in phase space — the instantaneous Lyapunov exponents — might become 
isotropic in the large-system limit. This suggests that the measured growth 
and decay rates are also closely associated with the subspaces spanned by 
the corresponding phase-space offset vectors {5}. Because the instantaneous 
Lyapunov exponents measure radial expansion and contraction, without any 
explicit rotational contribution, the instantaneous growth rates {\j(t) = Ajj} 
associated with every one of the phase-space directions contributing to a par- 
ticular vector are identical. The contributions in a fixed phase-space direction 
are the summed-up contributions from the entire set of Lyapunov vectors. 

At a particular phase-space point these contributions of the principal axes 
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of the dynamical matrix to the Lyapunov exponents are all proportional to 
cos 2 (#), where 9 is the angular difference between the principal axis and the 
corresponding instantaneous 5 vector. Because the Lyapunov exponents mea- 
sure growth or decay in the direction of S, the same logarithmic growth rate, 
(d/dt) In 5, applies to each component of the vector. 

In the full phase space the information dimension of the strange attractor is 
given by the number of Lyapunov exponents whose sum is zero: 

k 

]TA i = — ► D = k ; AD — AN + 2- k . 

3=1 

It has to be emphasized that k (as well as k^ introduced below) is not an 
integer. A linear interpolation between two successive values of the Lyapunov 
sum is implied, with k chosen such that the interpolated sum is precisely zero. 
In the Hamiltonian subspace the Lyapunov exponents contribute according to 
their projections into that space: 

E(cos 2 (^)A,(t)) = — D H = £<cos 2 (^-)> ; 

3=1 3=1 
4JV+2 k n 

AD n = £ (cos 2 (^)> = AN — 8 — £(cos 2 (^)> . 

3=kn 3=1 

The cos 2 (6?) form of the projection is required by the condition that chang- 
ing the sign of the offset vectors leaves the growth rate unchanged. This form 
follows from the quadratic form describing an infinitesimal phase-space hyper- 
ellipsoid centered on a moving trajectory point. Note also that this weighting 
satisfies the normalization of the projections: 

4JV+2 

]T cos 2 (^-) =AN-8 . 

3=1 

Partial sums give an effective number of exponents in the Hamiltonian projec- 
tion of the full phase space. Thus the analog of Kaplan and Yorke's conjecture 
for the Hamiltonian subspace is the effective number of exponents, ^(cos 2 (6)), 
at which the projected sum, XX-M^) cos2 (^))> vanishes. 

We explored the Kaplan- Yorke analog for the present problem, computing the 
delta vectors and their projections. The results are interesting. First, we noted 
that the summed-up local Lyapunov exponents are not strongly correlated 
with the directions in the subspace: 

EA(cos 2 W)^E(^)cos 2 W). 

This suggests that all the time-averaged projections of the various vectors into 
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Fig. 4. Time-averaged projections {{cos 2 (6))k} of the full-space {5k} into the Hamil- 
tonian portion of phase space for 4x4, 8x8, and 12 x 12 particles. The Nose-Hoover 
time r is 6. It is evident that in the large-system limit the influence of the boundary 
degrees of freedom disappears. 

the Hamiltonian subspace are similar, 



See Figs. 4 and 5. Fig. 4 demonstrates that the projection of the vectors 
becomes increasingly uniform as system size is increased, and is quite close to 
the average value, (AN — 8) /(AN + 2) expected for fully random projection 
directions. Fig. 5 compares the two estimates for the projected Lyapunov sums, 



The two sums vanish at nearly the same projected dimensionality, X)cos 2 (#), 
indicating that the correlation of the exponents with direction is small. 

We have used the analog of the Kaplan- Yorke formula to estimate the di- 
mensionality reduction in the Hamiltonian subspace, AZ}^, and show these 
results in Fig. 6. Just as in the full phase space, the loss of dimensionality 
varies smoothly with relaxation time and system size, and may be represented 



It is worth pointing out that a naive approach to dimensionality loss in the 
Hamiltonian subspace could be based on an orthonormalization of the Hamil- 



(cos 2 (#)) ~ (4iV - 8) /(AN + 2) . 



]TA(cos 2 (fl)> and £<A(t) cos 2 (fl)) . 



by 



AD n = 41.2(5) 



94(3) 10.1(5) 



12 




Fig. 5. Typical summed Lyapunov spectra for the full space and for two estimates 
of the Hamiltonian projection. The data shown here, for 64 particles and for a ther- 
mostat response time r = 5, show partial sums of the 258 Lyapunov exponents, 
52j< k Xj, and the two projections, J2j<k Aj(cos 2 (6>,,-)) ; ^2j< k (Xj(t) cos 2 (0j)) , as 
functions of k and X^<fc( cos2 (^j))- The two projections are actually different but 
indistinguishable on this scale. Kaplan and Yorke's form for the information dimen- 
sion is the linearly-interpolated value of k for which the linearly-interpolated sum 
vanishes. The analog, for the Hamiltonian portion of phase space, is the interpolated 
value of cos 2 (6 l ) at which the corresponding projected sum of exponents vanishes. 

tonian subspace only, propagating 5 vectors with the Hamiltonian equations of 
motion but using the thermostated equations of motion for the underlying ref- 
erence trajectory. This approach, which we explored years ago and which has 
recently been considered independently by Ken Aoki [private communications, 
2002], is equivalent to considering the Lyapunov spectrum for a Hamiltonian 
system subject to time-dependent forces {F(t)}, where the forces are due 
to the thermostated particles. Such equations of motion, being Hamiltonian, 
satisfy Liouville's Theorem [9,18], and guarantee that the corresponding Lya- 
punov spectrum is made up of pairs of positive and negative exponents, with 
zero sum. 



6 Conclusions 



The "0 4 " model shows conclusively that the nonequilibrium steady-state loss 
of phase-space dimensionality can easily exceed the dimensionality associated 
with the system boundaries. The present results also confirm that the de- 
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Fig. 6. Comparison of the dimensionality loss ADq-c in the Hamiltonian subspace 
to a fit, AD n = 41.2 - 94/\/iV - 10.1/r , to the data. 

viations from the large-system limit vary according to a power law, N~ p , 
l/4<p<l/2, roughly compatible with the inverse of the system size, 1/y/N 
in two dimensions. For r < 5 they have in addition deviations in the relaxation 
frequency oc (1/t). The quartic tethers are a particularly useful feature of the 
model, which make it possible to carry out the simulations without the need 
to take thermal expansion explicitly into account. 

The dimensionality loss AD is expected to be extensive as has been stressed 
in the introduction. For fixed thermostat temperatures T# and Tc the tem- 
perature gradient is determined by the system size, proportional to 1/ \/N in 
two dimensions, and AD becomes independent of N. This is indeed the case in 
the large-particle limit iV — > oo. Our simulations thus confirm the extensivity 
of the dimensionality reduction for stationary heat flow in the linear-response 
limit. Far from equilibrium weak deviations from this iV-dependence are found. 

Table 1 reveals that for short Nose- Hoover relaxation times the dimensionality 
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loss ADn in the purely Hamiltonian part of phase space may even slightly 
exceed the dimensionality loss AD in the full phase space. This is because 
AD-j-i is overestimated by one due to the lack of energy conservation in the 
Hamiltonian subspace once the thermostats are added. 

The projection of the phase-space offset vectors {6} into the Hamiltonian sub- 
space developed here shows that most of the dimensionality loss occurs in a 
part of the system which obeys purely conservative Hamiltonian equations 
of motion. Rapid rotation is responsible. This rotation nearly eliminates the 
correlations between phase-space contraction and direction. Evidently phase- 
space contraction is not only real, but relatively simple, and certainly must 
persist in the large-system thermodynamic limit. Thus the present results cor- 
roborate our interpretation of the second law of thermodynamics for nonequi- 
librium stationary states [8]. Such states occupy not just a reduced volume in 
phase space. They are restricted to a subspace of reduced dimensionality, with 
the dimensionality loss simply related to the rate of external thermodynamic 
entropy production. Useful models illustrating dimensionality loss for periodic 
color conducting or shear flows can now be developed as extensions of this 
idea. 
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